#include "../include/heat.hpp"

using namespace std;

int main()
{
    double h = 0.05;
    double k, r, nu = 1;
    r = 1;
    k = r * h * h / nu;
    int m; // nk(t) ih(x)
    VectorXd U;
    m = 1.0 / h;
    U = VectorXd::Zero(m + 1);
    for (int i = 0; i <= m; i++)
        U(i) = f1(i * h);

    MatrixXd A = MatrixXd::Zero(m + 1, m + 1);
    // set A
    for (int i = 0; i <= m - 1; i++)
    {
        A(i, i) = -2;
        A(i + 1, i) = 1;
        A(i, i + 1) = 1;
    }
    A(m, m) = -2;
    MatrixXd I = MatrixXd::Identity(m + 1, m + 1);
    double epsilon = 1e-6;

    int n = 0;
    while (n < 15)
    {
        VectorXd y1 = (I - 0.5 * r * A).lu().solve(r * A * U);
        U = U + y1;
        cout << U << endl;
        n++;
    }
}